Measuring spike train synchrony 
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Estimating the degree of synchrony or reliability between two or more spike trains is a frequent 
task in both experimental and computational neuroscience. In recent years, many different methods 
have been proposed that typically compare the timing of spikes on a certain time scale to be 
optimized by the analyst. Here, we propose the ISI-distance, a simple complementary approach 
that extracts information from the interspike intervals by evaluating the ratio of the instantaneous 
firing rates. The method is parameter free, time scale independent and easy to visualize as illustrated 
by an application to real neuronal spike trains obtained in vitro from rat slices. In a comparison 
with existing approaches on spike trains extracted from a simulated Hindemarsh-Rose network, the 
ISI-distance performs as well as the best time-scale-optimized measure based on spike timing. 
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INTRODUCTION 



The basic elements of neuronal communication are 
pulsed electric signals called action potentials or spikes. 
Under the assumption that both the shape of the spike 
and the background activity carry minimal information, 
neuronal responses are typically reduced to the much 
simpler form of a spike train, where the only informa- 
tion maintained is the timing of the single spikes. Mea- 
suring the overall degree of synchrony between different 
spike trains is an important tool in many different con- 
texts. It can be used to quantify the reliability of rc- 
spo nses upon repeated presentati ons of the same stimu- 
lus (jMainen and Seinowskil . Il995h . to address questions 
regarding the limitations o f neur onal coding (rate versus 
time coding) (jRieke et al. |, I1996D or to evaluate the in- 
formatio n transfer bet ween synaptically coupled neurons 
(cf.. e.g. jReved .l200i. 

A variety of different measures have been introduced 
to address the synchrony between spike trains. Most of 
these measures require considering a large number of tri- 
als, not just two. Some of them are based on the con- 
struction of a post-stimulus time histogram (PSTH), de- 
rived from multiple trials. PST H measures such as reli- 
ability, precision and sp arseness ( Mainen and Seinowskil . 
1995;' iBerrv et all 119971 ) rely on the analyst to define the 
so-called events, i.e., bursts of high firing rate. Other 
methods (i) quantify the occurrence of given spike pat- 
terns aiid_jnea^uretlreir_robustness ("attractor reliabil- 
ity" , iTiesinga et all l2002t ) , (ii) exploit the deviation of 
the spike train statist i cs from a Poiss onian distribution 
(|Brenner et all l200d : iTiesingal 12004 ) . or (iii) measure 
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the normalized var iance of poo l ed ex ponentially con- 
volved spike trains (jHunter et allfl998l ). 

The focus of this study lies on a group of measures 
that aim at a quantification of the degree of similar- 
ity or dissimilarity between as few as two spike trains. 
A very prominent example of such bivariate approaches 
are spike train distances that consider spike trains to 
be points in an abstract metric space and assign non- 
negative values quantifying the dissimilarity between a 
given pair of spike trains. Among these is the distance 
introduced in IVictor and Purpura! (1996), which evalu- 
ates the "cost" needed to transform one spike train into 
the other, using only ce rtain elemen t ary s teps. An- 
other metric proposed in Ivan Rossuml (|200lf) . measures 
the Euclidean distance between the two spike trains af- 
ter convolution of the spikes with an exponential func- 
tion. Other approaches quantify the cross correlation 
of spike trains a f ter exponential or Ga u ssian filtering 
(jHaas and White! . 120021 : ISchreiber et~aTl 120031) . or ex- 
ploit the exponentially weighted dis tance to the nearest 
neighbor (jHunter and Miltonl . 120031 ) . A common prop- 
erty of all these measures is the existence of one param- 
eter that sets the time scale for the analysis. This intro- 
duces elements of human fallibility and presumptions into 
the analysis and furthermore discourages direct compar- 
isons between different sets of results. Such a parameter 
does not exist for event synch r oniza tion, a method pro- 
posed in lOuian Quiroga et al.l ( 20021 ) that quantifies the 
number of quasi-simultaneous appearances, using a vari- 
able time scale that automatically adapts itself to the 
local spike rates. 

In this study a measure is proposed that, complemen- 
tary to the approaches mentioned above, uses the inter- 
spike interval (ISI) instead of the spike as the basic ele- 
ment of comparison. The ISI-distance quantifies the ratio 
of instantaneous rates and facilitates visualization of the 
relative timing of pairs of spike trains. Since no binning 
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is used, the measure has the maximum possible time res- 
olution (i.e., up to a single spike). Similar to event syn- 
chronization, it is both parameter-free and self-adaptive 
so that there is no need to fix a time scale beforehand. 
Thus the analyst is removed from the analysis, allowing 
for a more objective and broadly comparable measure of 
neuronal synchronization. In the first part of this study, 
the ISI-distance is illustrated using real in vitro data from 
cortical cells. 

Moreover, since a comparison between different ap- 
proaches was still missing, in the second part we tested 
the performance of several bivariate measures, includ- 
ing the ISI-disance, using clustered spike trains in a 
controlled setting. These were generated from a net- 
work of simulated Hindemarsh-Rose neurons with a pre- 
determined degree of coupling between pairs. In this sce- 
nario, different spike trains belonged to different clusters 
and the capability of the measures to detect the original 
clustering behavior could be quantified by two indices, 
which evaluate the correctness of the clusters and the 
separation between them, respectively. However, we do 
not claim to assess the absolute validity of our or any 
other measure, as no single number can represent the 
true code of the spike-generating mechanism under any 
circumstances. In a final step, to assess the similarity 
of the different approaches to measure spike train syn- 
chrony, we evaluated the degree of redundance between 
the different measures by means of a correlation analysis 

a 

The remainder of the paper is organized as follows: 
In the methods section, after a short description of the 
spike detection algorithm (Section III A[) . a more detailed 
description of the new method based on the ISI-distance 
is given (Section III Bp . It is illustrated using in vitro 
recordings from cortical cells in the entorhinal cortex 
of rats. The following section, III CI contains a short 
overview over the existing measures against which this 
new method is compared. The cluster analysis is de- 
scribed in Section [H Dl while the correlation analysis is 
described in section Hi El In section HlI Al the actual com- 
parison of the different methods on simulated time series 
taken from a network of Hindemarsh-Rose model- neurons 
is carried out. Conclusions are drawn in section HVl Fi- 
nally, both data sets (the in vitro recordings and the 
simulated Hindemarsh-Rose time series) are described in 
the appendix in sections lA II and lA 21 respectively. 



train can then be expressed as a series of <5 functions 



S(t) = £*(*-*i) 



(1) 



with ti, ...tM denoting the series of spike times and M 
being the number of spikes. 

In this study, for all the different measures the same 
spike detection algorithm is used. The threshold is chosen 
as the arithmetic average over the minimum and maxi- 
mum value of the action potential. 



B. The ISI-distance 

To obtain a time-resolved measure of the firing rate of 
the spike train {if}, in a first step the value of the current 
interspike interval is assigned to each time instant [32l | (see 
Figs. [13 top), 

x lsi {t) = min(if |£f > t) — max(tf |if < t) if <t < t x M 

(2) 

and accordingly for the second spike train {t?}. Now, 
in a second step the ratio between Xi S i and is taken 
(effectively, this is done just once after every new spike 
in either time series), and the final measure is thereby 
obtained after introducing a suitable normalization, 



7(f) 



Xi s i(t)/y isi (t) - 1 if x lsl (t) <= yui{t) 
-(yiai{t)/x iai (t) — l) else. 



( 3 ) 

The measure becomes zero in case of iso-frequent behav- 
ior, and approaches —1 and I, respectively, if the firing 
rate of the first (or second) train is infinitely high and 
the other infinitely low (see Figs. [HO bottom). 

Finally, in order to derive a measure of spike train dis- 
tance, there are two possible ways of averaging. In the 
time- weighted variant, the absolute ISI-distance is inte- 
grated over time, 



D I= I dt\I(t)\, 
't=o 



(4) 



whereas in the spike- weighted variant, the ISI-distance is 
evaluated only after every new spike in either time series, 



M 



II. 



METHODS 



A. Spike detection 

A prerequisite to any method is the extraction of the 
spike times from the time series by means of a standard 
spike detection algorithm. Typically, some sort of thresh- 
old criterion is employed, either for the time series itself 
or its derivative. Thereby the continuous time series is 
transformed into a discrete series of spikes. Each spike 



D s , 



(5) 



There are a number of meaningful extensions of this mea- 
sure. Omitting the absolute values yields a quantity that 
evaluates the relative firing rates of the two spike trains. 
Quantifying higher moments of the /(i)-distribution such 
as the standard deviation can provide additional infor- 
mation (in particular it allows to distinguish the cases of 
random jitter and systematic phase lag that could lead 
to the same ISI-distance). Also all of these variants can 
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Time [s] 

FIG. 1: (color online) First example of cortical cell recordings. 
In the middle traces the two recorded time series are shown. 
The detected spikes are marked in blue (input) and red (out- 
put), respectively. On top the ISI-values according to Eq. [2] 
are depicted, at the bottom the corresponding renormalized ISI- 
distance (cf. Eq. [4}. Here colors mark the times where the 
respective spike train is slower. For this pair of spike trains an 
ISI-distance Di = 0.06 is obtained. 



be implemented using a moving- window analysis. Due 
to the self-adaptation of this measure to the time scale 
of the spikes, reasonable results can be obtained also for 
rather short spike trains. Finally, the sensitivity of the 
measure can be extended to longer time scales and a more 
coarse-grained evaluation by averaging Xi S i (t) and yi S i (t) 
over neighboring ISIs. An example of a possible appli- 
cation is the quantification of (dis-)similarities between 
bursts in time series. 

In Fig. []] the ISI-distance is applied to two exemplary 
input-output spike trains of 10 s duration (for a descrip- 
tion of the data see Appendix I A In the first seconds, 
the spike trains are 1 : 1 synchronized and this is re- 
flected by an ISI-distance I(t) 0. Nevertheless, small 
deviations can be visualized that are hard to catch from 
a visual inspection of the spike trains themselves. These 
deviations are more pronounced in the second half of the 
recording where the output no longer follows the input 
but rather slows down (as reflected by predominantly 
negative values marked in red) and a spike doublet occurs 
(as indicated by the short excursion to positive values 
marked in blue) . The example shown in Fig. [2] reveals 
that certain patterns in the ISI-distance appear repeat- 
edly. The output exhibits several spike doublets, some of 
which arc followed by a miss (reflected by the negative 
values marked in red). Finally, a more irregular behavior 
is shown in Fig. [3] where it is again clear that the ISI- 
distance allows tracing the relative firing rate behavior 
in a simple way [33| . 
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FIG. 2: (color online) Second example of cortical cell recording 
shown in the same way as in Fig. [T] In this case the ISI-distance 
attains the value Di — 0.11. 
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FIG. 3: (color online) Third example of cortical cell recording 
shown in the same way as in Fig. [TJ The ISI-distance for this 
example is Di = 0.34. 



C. Existing measures of spike train distance 



In this study the ISI-distance will be compared 
against five existing measures of spike train (dis- 
similarity. These comprise the s pike t rain metrics in- 
troduced in IVictor and Purpura! (199(f), as well as in 



Ivan Rossuml (l200ll). a correlation measure proposed in 
ISchreiber et al.l (120031). another distance measure intro- 



duced in Hunter and Miltonl ( 20031). and event synchro- 
nizati on, a method introduced in iQuian Quiroga et alj 

i|2oo2h [a3- 

In order to compare the various measures, we turned 
each of them into a suitably-normalized measure of dis- 
similarity, as this is more akin to the concept of distance. 
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1. Victor- Purpura spike train metric 

The spike tra in me tric Dy introduced in 
IVictor and Purpura! ( 19961 ) defines the distance be- 
tween two spike trains in terms of the minimum cost of 
transforming one spike train into the other by using just 
three basic operations: spike insertion, spike deletion 
and spike movement. While the cost of insertion or 
deletion of a spike is set to one, the cost cy of moving 
a spike is the only parameter of the method setting the 
time scale of the analysis. For small cy, the distance 
basically equals the difference in spike number, whereas 
for high cy, the distance approaches the number of 
non-coincident spikes, since instead of shifting spikes 
it becomes more favorable to delete all non-coincident 
spikes of the one time series and to insert all non- 
coincident spikes of the other. Thus, by increasing the 
cost cy , the distance is transformed from a rate distance 
to a timing distance. 



2. Van Rossum spike train metric 

A second spike train metric was introduced in 
Ivan Rossu m (2001). In this method, each spike is con- 
volved with an exponential function e~^ t ' ti ^ TR (t > ti), 
where ti is the spike time. From the convolved wave- 
forms f(t) and g(t), the van Rossum distance Dr can be 
calculated as 



D R {r R ) = - [°°lf(t) 

T R Jo 



g(t)] 2 dt 



(6) 



Since the post-synaptic currents triggered by the single 
spikes approximate exponentials, the van Rossum dis- 
tance estimates the difference in the effect of the two 
trains on the respective synapses. In this method, the 
time constant tr of the exponential is the parameter set- 
ting the time scale. It is the inverse of Victor and Pur- 
pura's cost parameter: tr — l/cy. 



3. Schreiber et al. similarity measure 

The correlation-b as ed ap proach was first described 
in lHaas anefw hite (2002), and later detailed in 
ISchreiber et al.l (|2003l) . In this approach, each spike train 
is convolved with a filter o f a certain wi dth (exponential 
in lHaas and White! (|2002h . Gaussian in I Schreiber et al.l 
(|2003D ) to form Si before cross correlation and normal- 
ization. 



Ss(vs) 



* 3 



(7) 



Haas and White allowed a minimal phase lag in the cross 
correlation (and thus another parameter to adjust), while 
Schreiber et al. allowed none. Here, the approach by 
Schreiber et al. is used. The width of the convolving 



filter as sets the time scale of interaction between the 
two spike trains. The inversion Dg = 1 — S$ yields a 
normalized measure of spike train dissimilarity that can 
be compared with the other. A cl ustering analy s is bas ed 
on this measure was performed in lFellous et al.l (|2004l ). 



4- Hunter-Milton similarity measure 



Thi s approach, introduced in iHunter and Milton! 
(2003), starts by identifying in the spike train y the near- 
est spike tf.(f. to the spike occurring at time if in the 
spike train x. The degree of coincidence between these 
spikes is quantified by r xy = exp (— \t® — ^^I/Tffj an d 
the overall similarity measure Sh is thereby determined 
as the symmetrized average of r xy over the entire series, 



Sh = 



(8) 



Also in this method, there is a free time scale that can be 
set by fixing the parameter th- For two identical spike 
trains r xy = r yx = 1. Accordingly, a measure of spike 
train dissimilarity can be obtained as Dh = 1 — Sh 



5. Event synchronization 

The last measure i s a variant of the ev e nt sy nchro- 
nization pr oposed in Quian Quiroga et al. (l2002f) . an d 
later used inlHahnloser et all (|2002l ). lKreuz et al.l (j2004h , 
and iKreuz et al. I (120071) . This measure quantifies the 
overall level of synchronicity from the number of quasi- 
simultaneous appearances of spikes. However, in con- 
trast to the measures introduced above, this method is 
scale-free, since the maximum time lag Ty up to which 
two spikes tf and $ are considered to be synchronous is 



adapted to the local spike rates according to 

mi,:{/;. ; /;•/; /; '? , I 2. (9) 

Then the function c(x\y) is introduced to count the num- 
ber of times a spike appears in x shortly after a spike 
appears in y, 



M x My 

c{x\v) = ^2Y1 J v> 
i=i j=i 



(10) 



where 



Ji 



if < tf - 
if t x = t y 
else. 



t v < t- 



(11) 



With c{y\x) defined accordingly, we can write the event 
synchronization as 



Q 



c(y\x) +c(x\y) 

y/M x M y 



(12) 



5 



Again, a measure of spike train distance can be denned as 
Dq = Q— 1. With the above normalization, < Dq < 1, 
with Dq = if and only if all spikes of the signals are 
synchronous [35[. 



D. Assessing clustering quality 

One important application for measures of spike train 
synchrony is the identification of interspike correlations 
and the reconstruction of clustering patterns. In order to 
test the above measures, we generated 29 spike trains by 
simulating a network of Hindmarsh-Rose neurons (see ap- 
pendix |X2|. From the network architecture and the pat- 
tern coding, we organized the 29 spike trains into three 
principal clusters with 13 members (clusters 1 and 2) and 
3 members (cluster S), respectively. We first validated 
the different measures and then quantified their perfor- 
mance in reproducing the cluster structure by means of 
two indices. For the four measures Dy, Dr, D$ and Dh 
that depend on a parameter setting the time scale, we 
varied the respective parameter over several logarithmic 
decades, with four equidistant values within each decade, 
i.e., cy = io _a + - 25h . In each case, we adapted the pa- 
rameter range via a and b to cover the relevant extreme 
cases. 

After applying a given similarity measure to all pos- 
sible pairs of spike trains, we generated a hierarchical 
cluster tree (dendrogram) by applying the single linkage 
algorithm provided by Matlab to the resulting pairwise 
distance matrices. An exemplary dendrogram obtained 
from the 29 Hindemarsh-Rose time series using the event 
distance Dq is shown in Fig. 2] Three principal clusters 
are clearly distinguished. The dendrogram is constructed 
as follows: First, the closest pair Si, Sj of sequences is 
identified and thereby linked by a n-shaped line, where 
the height of the connection measures the mutual dis- 
tance d(Si, Sj). These two time series are merged into a 
single element C a , and the next closest pair of elements is 
then identified and connected. The procedure is repeated 
iteratively until a single cluster remains. The implemen- 
tation of the method requires introducing the distance 
between a pair of clusters C a , Cp. In the single linkage 
algorithm, it is defined as the minimum over all the dis- 
tances between pairs of sequences in the two clusters, i.e., 
d(C a ,Cp) = mm{d(S k ,S m )}, S k G C a , S m G Cp. 

In order to quantify the success in reproducing this 
clustering, we comput e d the entropy of the confus i on ma - 
trix N af3 ) (|Abramsonl Il963t IVictor and Purpural 1 19961) . 
The entry N a p is defined as the number of times Sp is the 
clo sest cluster to a spike train belonging to S a . Follow- 
ing IVictor and Purpural the distance between the spike 
train Si and the cluster C a is defined as (d(Si, Sj)) a , 
where (-)j denotes the average over all spike trains in the 
cluster C a . Note that this distance is different from the 
one used for the cluster identification, however, we veri- 
fied that results proved to be robust against variations of 
the distance used. For a perfect clustering N is diagonal, 
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FIG. 4: (color online) Example of a hierarchical cluster tree 
obtained from 29 Hindemarsh-Rose time series employing the 
event distance Dq. The three principal clusters C\, C2 and Cs 
are distinguished by different colors. The merging of the first 
two time series Si and Sj to cluster C a is highlighted by very 
thick blue lines, the consecutive merging of this cluster with Sk 
by thick blue lines. Finally, black lines mark the separation of 
the three principal clusters used for the definition of the cluster- 
separation F. The clustering performance values for this example 
are H = 1 and F = 0.57. 



whereas each misclassification yields a non-diagonal ele- 
ments. The relative amount of misclassifications is finally 
quantified by the entropy 



tt \ ~* 1 Pap 

c= 5 ^ 



(13) 



where p a p = N a p/N ta t, P a = J2bPab, and Qp = J2bPbfS- 
This entropy value is then normalized to the maximum 
entropy obtained for a correct classification [36| . 



H — Hc/H„ 



(14) 



Although the clustering entropy H evaluates the correct- 
ness of the hierarchical tree, it does not quantify the sepa- 
ration of the three principal clusters in those cases where 
the expected clustering is obtained. As can be seen in 
Fig. |H the cluster separation is given by the lengths L\, 
L2, Ls, and Lm of their upper branches. A convenient 
way of quantifying the cluster separation with a single 
indicator is by introducing the index 



F 



Li+ L 2 + L s + L 



M 



3L A 



(15) 



where the branch lengths are normalized to the difference 
La between the overall maximum and the overall mini- 
mum distance thus guaranteeing that the F-values range 
in the interval [0, 1] . Given a correct clustering with three 
principal clusters, all quantities needed for the calcula- 
tion of F can be extracted from the output matrix of the 
Matlab single linkage algorithm, otherwise the clustering 
separation is set to F = 0. 



E. Correlations between the different measures 

In order to investigate to which extent the different 
measures of spike train distance carry independent and 
non-redundant information, we performed a correlation 
analysis. First, for each of the four measures Dy, Dr, 
Ds, and Dh that depend on a parameter setting the 
time scale, we identified the parameter value yielding 
the best clustering results. Then, we applied each of 
the six measures (Dj, Dy, Dr, Ds, Dh, and Dq) to 
the different pairs of sequences, obtaining six sets of 
29 x 28/2 = 406 different values. In order to guar- 
antee maximal homogeneity, the various measures were 
all scaled to the [0, 1] range (this means that the two 
unnormalized measures Dy and Dr were divided by 
their maximal values). Moreover, following a custom- 
ary approach to better fit a normal distribution, each 
data set was arcsin-t r ansfq rmed using x' = arcsin(y / x) 
([Daniels and Kendalll . fl947T) . Fi nally, we determined the 
Pearson correlation coefficients ( Devore and Peckl 120051 ) 
and from the pairwise distances (1-correlation) among 
the different measures, we obtained a hierarchical cluster 
tree (using again the single-linkage algorithm) . 



III. 



RESULTS 



Comparison of measures using simulated 
Hindemarsh-Rose time series 



In order to compare the various dissimilarity mea- 
sures, we have analyzed numerically generated spike 
trains, since their properties are much more controllable. 
More precisely, we refer to Hindemarsh-Rose time series 
whose clustering properties are known beforehand (see 
appendix I A 21 for a description of the underlying model). 
Two instances of the ISI-distance are shown in Figs. [5] 
and [5] where the signals emitted by two neurons belong- 
ing to the same and to different clusters, respectively, are 
compared. In the first example, deviations from zero of 
the ISI-distance arc confined to short time scales (they 
are mostly due to small phase shifts that accompany large 
changes of the firing rate). In the second example, long- 
lasting differences are detected which also exhibit a sort 
of oscillation. 

The distance matrix obtained from the application of 
the ISI-distance Dj to all 406 combinations of the 29 
spike trains is shown in Fig. [7] Patterns can be clearly 
recognized, since the neurons have been ordered accord- 
ing to to their a priori affiliation known from the model 
setup. Smallest distances are obtained for pairs of spike 
trains belonging to the same cluster starting from those 
within S. At the other extreme, the largest distances are 
found for spike trains belonging to the clusters 1 and 2. 

From this distance matrix, we generated the cluster 
dendrogram shown in Fig. [5J where the three principal 
clusters are clearly visible. The absence of misclassifica- 
tions implies that the clustering entropy (computed ac- 
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FIG. 5: (color online) Two time series from two neurons coding 
both for the same pattern B (middle). The detected spikes are 
marked in blue and red, respectively. On top the I S Is, at the 
bottom the corresponding renormalized ISI-distance. Here colors 
mark the times where the respective spike train is slower. For 
this combination an ISI-distance Di — 0.019 is obtained. 




Data points 



FIG. 6: (color online) Same as Fig. [5]but this time for two time 
series from two neurons coding for different patterns. In this 
case the ISI-distance Dj = 0.032 is much higher. 



cording to Eqs. (|13ll4j) ) is H = 1. On the other hand, 
the cluster separation determined from Eq. I|15p is for 
this case F — 0.66. Similar results have been obtained 
for the other parameter-free measure, the event distance 
Dq (the corresponding dendrogram being shown in Fig. 
2]). Also in this case, H = 1, while the smaller value of 
F (0.60) suggests a slightly lower clustering quality. 

We then investigated the performance of the Victor- 
Purpura and van Rossum spike train distances, as well 
as of the Schreiber et al. and Hunter-Milton dissimilar- 
ities for different values of the free parameter in order 
to select the proper time scale. In Fig. |H] the Victor- 
Purpura spike train distance Dy is plotted against the 
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FIG. 7: (color online) Distance matrix for 29 time series from the 
Hindemarsh-Rose model. Results are obtained by using the ISI- 
distance Di . Neurons are labelled by '1', '2' and 'S', respectively, 
depending on their affiliation to pattern 1, 2 or both ('shared'). 

0.3 h 
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Sa Sb Sc 2a 2e 2m 2b 2j 2i 2k 2d 2h 21 2f 2g 2c 1a 1c 1d 1i 1j 1b 1h 1m 1f 11 1g 1k 1e 
Neuron 

FIG. 8: (color online) Clustering of 29 time series from the 
Hindemarsh-Rose model using the ISI-distance. The existence 
of three clearly separated clusters is evident. The cluster perfor- 
mance values are H — 1 and F — 0.66. 

cost parameter cy for all pairs of spike trains in a group 
containing three members in each of the three principal 
clusters. The relative order of the six different combina- 
tions of cluster affiliations (1 — 1, 1 — 2, 1 — 5, 2 — 2, 
2 — 5, 5 — 5) depends on the cost parameter. At small 
cy values, the Victor-Purpura distance measures the dif- 
ference in spike counts and this number seems not to be 
closely related to the type of underlying cluster (see in 
particular the spread of the Dy-values corresponding to 
the 1 — 2 combination). At high cy values, distances are 
larger and also very mixed. It is only at intermediate 
values that a clear separation of the six different cluster 
combinations can be observed. High values are obtained 
for all the 1 — 2 combinations while the intra-cluster 5—5 
distance is the smallest one. 

We find similar results for the other parameter- 
dependent measures. In all cases, there exists an inter- 
mediate parameter range where meaningful results can 
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FIG. 9: (color online) Dependence of the Victor-Purpura spike 
train distance Dv on the cost parameter cy . For each combi- 
nation of cluster affiliations a different color is used. The dotted 
line marks the cost value for which the best clustering separation 
is obtained. On the right hand side of the dashed line the values 
obtained for the ISI-distance without any optimization are de- 
picted. Note that due to different normalizations and scalings, 
the measures appear on different y-axes. 



be obtained (i.e., the cluster entropy is equal to 1), while 
for higher and lower values no clear clustering can be 
recognized. The role of the free parameter is better seen 
by determining the clustering performance for different 
values of the parameter itself. In Fig. [TU] the two indices 
H and F are plotted for the Victor-Purpura distance Dy 
versus the cost parameter cy . We see that a perfect clus- 
tering is found only inside an intermediate range, where 
H = 1; the smaller H values found outside this interval 
reflect the presence of misclassifications in the clustering 
tree. Inside the interval where the right classification is 
obtained, we computed the clustering separation F . This 
index attains its peak value F = 0.67 when cy = 0.01, 
which we thus identify as the optimal value of the time 
scale for the separation of the different clusters. This re- 
sult is consistent with the visual impression when looking 
at Fig. (the vertical line corresponds to the optimal cy 
value). A similar scenario is obtained also for the other 
measures that depend on a time scale. In each case there 
exists an intermediate range where H = 1. The broad- 
est range is found for the Victor-Purpura and the van 
Rossum distance. 

The performance of the different measures are com- 
pared in Fig. 1 1 X I where the maximum value of the clus- 
ter separation is shown for each spike train distance (for 
the parameter-free ISI-distance and the event distance no 
optimization is required). The best results are found for 
the Victor-Purpura distance Dy, but the ISI-distance Dj 
and the van Rossum distance Dr perform almost equally 
well. At the other extreme, the poorest cluster separation 
is obtained for the Hunter-Milton dissimilarity Dh ■ 
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FIG. 10: For the Victor-Purpura spike train distance the cluster- 
ing entropy H (solid line) and the cluster separation F (dashed 
line) in dependence of the parameter ey that sets the time scale. 
Remember that F is only calculated for those cases where a cor- 
rect clustering (as reflected by H — 1) is obtained. The optimum 
performance F = 0.67 is marked by a large asterisk at cy = 0.01. 
On the right hand side of the dashed line the respective values 
for the ISI-distance (H — 1 and F = 0.66), which required no 
parameter optimization, are depicted. 
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FIG. 12: (color online) Correlation coefficients between the six 
measures of spike train dissimilarity. 



FIG. 13: Clustering for the six measures of spike train dissimi- 
larity. 



B. Correlations between the different measures 

In order to assess the degree of redundancy among the 
different measures of spike train dissimilarity, a correla- 
tion and cluster analysis has been performed on all 406 
bivariate combinations of measure results by following 
the approach discussed in the previous section. Since 
the number of independent observations can hardly be 
estimated, this is only a relative examination. For this 
reason no values of significance are given. 

As we see from the high overall level of correlation 
shown in Fig. I12[ all measures seem to carry similar in- 
formation. The minimum correlation coefficient 0.91 is 
obtained between the ISI-distance and the Hunter-Milton 
dissimilarity, while the spike train distances by Victor- 
Purpura and van Rossum appeared to be the most corre- 
lated measures. From the corresponding cluster tree (cf. 



Fig. [T3"|) . it becomes evident that the ISI-distance is the 
most independent measure, whereas the other measures 
belong to one big cluster. This reflects the fact that the 
ISI-distance is derived from interspike intervals, while the 
other measures are based on spike times. 



IV. DISCUSSION 

In this study we propose a simple method for measur- 
ing the (dis-)similarity of two spike trains. As an esti- 
mator based on the relative sizes of interspike intervals, 
the ISI-distance is complementary to all spike-based mea- 
sures of synchrony that quantify the simultaneous occur- 
rences of spikes via some sort of coincidence detection. 
As illustrated by an application to in vitro recordings 
of cortical cells, the measure serves also as an excellent 
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means to visualize the occurrence of spiking pattern in 
the respective spike trains. Finally, this approach repre- 
sents a natural starting point towards a more complete 
characterization of neuronal activity, in so far as higher 
moments (such as the standard deviation) of the I(t) dis- 
tribution can be measured and compared. 

In order to judge the relative merit of the different 
methods, we compared the methods by evaluating each 
of them on spike trains extracted from a network of sim- 
ulated Hindemarsh-Rose neurons. We assessed the abil- 
ity of the different measures to reproduce the original 
clustering (established by a priori adjusting the synap- 
tic couplings in the model) by means of two indices. In 
this comparison, no measure fails in reproducing the ex- 
pected clustering; however, we found subtle differences in 
the degree of separation among the three clusters. The 
ISI-distance performed as well as the best spike-based 
measure (the Victor-Purpura distance Dy) with the dis- 
tinct advantage of not requiring the optimization of any 
parameter. In fact, the ISI-distance, like event synchro- 
nization, is self-adaptive in that it automatically identi- 
fies the proper time scale. In particular, this holds true 
for changes of firing rate within the same spike trains. 
Whereas the ISI-distance can also be applied to spike 
trains that include different time scales (i.e., regular spik- 
ing and bursting) the other measures would misrepresent 
either behaviour depending on the parameter chosen. 

Finally, we implemented a correlation analysis in order 
to evaluate the degree of independence among the six dif- 
ferent measures of spike train similarity. Since the overall 
level of correlation is quite high, all measures apparently 
access similar information. The subsequent cluster anal- 
ysis reveals that the ISI-distance is the most independent 
approach. This is not surprising since the ISI-distance is 
the only measure that can be regarded as a measure of 
rate coding (since it is built on instantaneous firing rate 
estimates) whereas all other measures (that are based 
on spike timing) can be interpreted as measures of time 
coding. 

Some general remarks concerning the application of 
measures of spike train (dis-)similarity: First, although 
the focus of this study is on the estimation of similar- 
ity between just two spike trains, all methods can also 
be used within a multivariate context. For example, in 
order to assess the reproducibility of neural spike train 
responses to an identical stimulus across many differ- 
ent presentations (trials), reliability can be defined as 
th e average over all pairw ise correlation-val ues (as done 
in lHaas and White! (l2002l) . ISchreiber etal] (|2003l) . and 
iHunter and Miltonl (|2003l )1. O n the other hand, it should 
be noted that the ISI-distance, like all the other mea- 
sures, can be fooled by phase lags, no matter whether 
these are caused by internal delay loops in one spike train 
or by a common driver that forces the two time series 
with different delays. Thus any such phase lags should 
be removed by suitably shifting the time series before 
computing the measures. 

The comparisons carried out in this study constitute an 



important step towards the identification of the most ap- 
propriate spike train (dis-)similarity for an application to 
real data. From the observed performance values and the 
results of the correlation analysis, we conclude that the 
ISI-distance together with either one of the spike based 
measures, such as the Victor-Purpura or the van Rossum 
spike train distance, or event synchronization is the most 
appropriate choice. However, it should also be noted that 
other possible criteria for the selection of an appropri- 
ate measure, such as computational cost or robustness 
against noise, remain to be evaluated in future studies. 



APPENDIX A: DATA 



1. Application to in vitro recordings of cortical 
cells 



The ISI-distance is illustrated using in vitro whole- 
cell recordings taken from cortical cells from the layer 
2 medial cntorhinal cortex of young Long-Evans rats. In 
these experiments (conducted as approved by the UCSD 
IACUC) cortical cells were selected from a slice prepara- 
tion (400 micron) by their superficial position, as well as 
particular characteristics of thei r electrophysiological re - 
sponses to long current steps (cf. lHaas and Whit el l2002D . 
Intracellular signals were amplified, low pass filtered, and 
digitized at 10 kHz via software created in Lab View (Na- 
tional Instruments). Inputs were delivered as synap- 
tic conductances th rough a linux-based dynamic clamp 
( Dorval et al.l . l200lh . Inputs were comprised of synaptic 
inputs added to an underlying DC depolarization. The 
amplitude of the DC depolarization was tailored for each 
cell to elicit a spike rate of 5 — 10 Hz. Synaptic inputs 



were of the form Isyn = G syn S(y m — V syn ) where S fol- 
lows the differential expression dS/dt = a(l — S) — j3S; 
a = 500/ms, = 250/ms, V syn — OmV. G syn was tai- 
lored for each cell to be peri-threshold and elicit a spike 
with probability close to 50%. Synaptic events were de- 
livered for 10 seconds at either regular intervals or chaotic 
intervals; the latter intervals were taken from a set of 
five 60 second recordings of SC spike times in response 
to steady DC depolarization alone. 



2. Hindemarsh-Rose simulations 

The spike trains have been generated using time se- 
ries extracted from a larger network of Hindemars h-Rose 
(HR) model-neurons (|Hindmarsh and Rosel Il984h in the 
chaotic regime. This network was originally designed to 
analyze semantic memory representations using feature- 
based models; details of the network architecture and 
th e implementa t ion of the feature coding can be found 
in lMorelli et~aTl (|2006t) . 

In short, the state the neuron i is determined by three 
first-order differential equations describing the evolution 
of the membrane potential Xi, the recovery variable V;, 
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and a slow adaptation current Z^ 

X t = Y t - Xf + 3Xf -Zi + Ii + Oi(t) - ft(t)(Al) 
% = \-hXf-Yi (A2) 
Z % = 0.006[4pQ - 1.6) - Zi), (A3) 



where 



F(M-l) 



(A4) 



kept constant, the 29 time series to be analyzed were ex- 
tracted. According to their coding properties regarding 
the retrieval of two distinguished memory patterns, the 
29 time series belonged to three principal clusters, 13 of 
the corresponding neurons coded for pattern 1 only, 13 
coded for pattern 2 only and 3 coded for both patterns 
(shared). The respective time series were labelled by T', 
'2' and 'S' followed by an index letter. The numerical 
integration was done by using a fixed-step fourth-order 
Runge-Kutta method. The integration step-size was cho- 
sen equal to 0.05 ms of real time. The length of the time 
series analyzed was 32768 data points. 



and 



A(*) = ^E A * ) (*)- 

fe=l 



(A5) 



The network consisted of N = 128 HR neurons be- 
longing to M = 16 different modules with F = 8 neu- 
rons each. In a learning stage, input memory patterns 
were stored by updating the synaptic connection weights 
Wij between different neurons using a Hebbian mecha- 
nism based on the activity variables Aj. During the re- 
trieval stage in which the learned connection weights were 
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